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Abstract 

Universal role of the nonlinear one-third subharmonic resonance 
mechanism in generation of the strong fluctuations in such complex 
natural dynamical systems as global climate and global solar activity is 
discussed using wavelet regression detrended data. Role of the oceanic 
Rossby waves in the year-scale global temperature fluctuations and the 
nonlinear resonance contribution to the El Nino phenomenon have 
been discussed in detail. The large fluctuations of the reconstructed 
temperature on the millennial time-scales (Antarctic ice cores data for 
the past 400,000 years) are also shown to be dominated by the one- 
third subharmonic resonance, presumably related to Earth precession 
effect on the energy that the intertropical regions receive from the 
Sun. Effects of Galactic turbulence on the temperature fluctuations 
are discussed in this content. It is also shown that the one-third 
subharmonic resonance can be considered as a background for the 11- 
years solar cycle, and again the global (solar) rotation and chaotic 
propagating waves play significant role in this phenomenon. Finally, 
a multidecadal chaotic coherence between the detrended solar activity 
and global temperature has been briefly discussed. 
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1 Introduction 



Since the global climatology is in an earlier state of data accumulation, get- 
ting knowledge from the data is not an easy and straight forward process. 
However, there is only one way enabling to understand physics of the phe- 
nomenon: accurate and punctual analysis of data. Moreover, the global 
climate, as we know it now, is a nonlinear phenomenon. This makes the 
process even more difficult and painful than usual. It is normal for this stage 
that the models are at best crude while limitations are great and (that even 
more dangerous) unknown. 

Already the simplest energy balance model for the space averaged surface 
temperature T has a form of a nonlinear equation: 

f = ^{^-*{T))-eaT*}) (1)) 

where C is the heat capacity. The first term in the right hand side of the 
balance equation represents incoming solar energy with q as the solar constant 
and a as the albedo, and second term represents outgoing infrared energy 
with e and a as the emissivity factor and the Stefan constant respectively. 

A more complex nonlinear model describing an interaction between sur- 
face energy balance and mass balance of the cryosphere has been suggested 
in Ref. [1] (in dimensionless excess variables): 

A = ciC, + c 2 s - s 3 - s 2 C + F sm(ut) (3) 

where s the sea-ice extent, ( = T — s and T is the mean ocean surface 
temperature, c\ and c 2 are certain constants, and F sm(ut) stands for solar 
forcing. It is shown in Ref. [2] that the system Eqs. (2-3) can be considered 
as a the perturbations (for 1 ^> s) of a reference system 

dx 

Tt =y (4) 

d V 3 (r\ 

— = cx — x (5) 
dt y ' 

which is a Hamiltonian system for the conservative Duffing oscillator. 
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The nonlinear Duffing oscillator appears also in the conceptual model of 
the tropical Pacific oscillations (for positive upper ocean temperature anoma- 
lies T driving shallow thermocline depth anomalies h) [3], [3] 



The periodic term in this equation represents the solar forcing (77 is an am- 
bient noise, and 7, h and F are certain geophysical constants, see for more 
details below). 

The appearance of the nonlinear Duffing oscillator in the models of the 
global climate processes for different time-space scales: from Paleoclimate 
Eqs. (2-5) to the decadal and inter-year oscillations Eqs. (6-7), shows that 
this nonlinear oscillator can be considered as a simple prototype of nonlin- 
ear climate systems (it will be also shown below how the nonlinear Duffing 
oscillator appears naturally in the solar dynamo models Eq. (17)). Basic 
symmetries of the systems can play a crucial role in this phenomenon (see 
below). For the non- linear systems far from equilibrium the oscillatory be- 
havior can appear near bifurcation points (the dynamics generated by the 
Duffing oscillator can exhibit both deterministic and chaotic behavior [5]- [7J, 
depending on the parameters range). 

One of the specific non-linear properties of the Duffing oscillator is the so- 
called one-third subharmonic resonance [8j. Let us imagine a forced excitable 
system with a large amount of loosely coupled degrees of freedom schemat- 
ically represented by Duffing oscillators with a wide range of the natural 
frequencies ujq (it is well known [9] that oscillations with a wide range of 
frequencies are supported by ocean and atmosphere, cf. also Ref. [TO]): 



where x denotes the temporal derivative of x, (3 is the strength of nonlinearity, 
and F and u are characteristic of a driving force. It is known (see for instance 
Ref. [8]) that when u> w 3a;o and (3 <C 1 the equation (8) has a resonant 
solution 




(6) 





x + UqX + 7X + f3x 3 = F sin cut 



(8) 
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where the amplitude a and the phase tp are certain constants. This is so- 
called one-third subharmonic resonance with the driving frequency u. For 
the considered system of the oscillators an effect of synchronization can take 
place and, as a consequence of this synchronization, the characteristic peaks 
in the spectra of partial oscillations coincide [TT| . It can be useful to note, 
for the global climate modeling, that the odd-term subharmonic resonance is 
a consequence of the reflection symmetry of the natural nonlinear oscillators 
(invariance to the transformation x — » — x, cf. also Ref. [12]). 

This non-linear phenomenon can be observed directly in the real data. 
Let us start from the well known (in relation to the global warming) graph 
representing the monthly global temperature data: figure 1 (the instrumen- 
tal monthly data are available at http : / / lwf . ncdc . noaa . gov /oa/ climate / 1 re- 
search/anomalies/index, html, see also Ref. [13]). One can see that the tem- 
perature is strongly fluctuating. Actually, the fluctuations are of the same 
order as the trend itself (see figures 1 and 2). While the nature of the trend is 
widely discussed the nature of these strong fluctuations is still quite obscure. 
As we will see further these strong fluctuations have an intrinsic nonlinear 
nature. This is the one-third subharmonic resonance to annual solar forcing. 
Cyclic forcing, due to astronomical modulations of the solar input, right- 
fully plays a central role in the long-term climate models. Paradoxically, it 
is a very non-trivial task to find imprints of this forcing in the long-term 
climate data. It will be shown in present paper that just unusual properties 
of nonlinear response are the main source of this problem. 

There are many well known reasons for asymmetry in response of the 
North and South Hemispheres to solar forcing: dominance of water in the 
Southern Hemisphere against dominance of land in the Northern one, topo- 
graphical imbalance of land (continents) and oceans in the Northern Hemi- 
sphere due to continental configuration, seasonality and vegetation changes 
are much more pronounce on land than on ocean surface, and anthropogeni- 
cally induced asymmetry of the last century. This asymmetry results in 
annual asymmetry of global heat budget and, in particular, in annual fluc- 
tuations of the global temperature. Nonlinear responses are expected as a 
result of this asymmetry. 

The nonlinear one-third subharmonic resonance produces strong fluctua- 
tions of the global temperature also for much larger paleoclimate time scales: 
figures 3 and 4 (the reconstructed data are available at |http:/ /www.ncdc 
noaa.gov/paleo/metadata/noaa-icecore-6076.html, see also Ref . [H]). Again, 
while the nature of the trend is widely discussed (in relation to the glacia- 
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Figure 1: The monthly global temperature data (dashed line) for the pe- 
riod 1880-2009. The solid curve (trend) corresponds to a wavelet (symmlet) 
regression of the data. 

tion cycles) the nature of these strong fluctuations is still quite obscure. In 
this case the one-third subharmonic resonance, generating the fluctuations, 
is presumably related to Earth precession effect on the energy that the in- 
tertropical regions receive from the Sun. 

We will also see that in the Sun itself Nature uses the nonlinear one-third 
subharmonic resonance to amplify the hydromagnetic dynamo effects to the 
observed 11-year periodic solar activity. And again, a global rotation and 
nonlinear waves are the main components of the nonlinear mechanism. 

Such universality of this non-linear mechanism makes it a very interest- 
ing subject for a thoughtful investigation. The problem has also a technical 
aspect: detrending is a difficult task for such strong fluctuations. Most of the 
regression methods are linear in responses. At the nonlinear nonparametric 
wavelet regression one chooses a relatively small number of wavelet coeffi- 
cients to represent the underlying regression function. A threshold method 
is used to keep or kill the wavelet coefficients. At the wavelet regression the 
demands to smoothness of the function being estimated are relaxed consider- 
ably in comparison to the traditional methods. These advantages make the 
non-linear wavelet regression method an adequate tool for detrending the 
strongly fluctuating natural data. 
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Figure 2: The wavelet regression detrended fluctuations from the data shown 
in Fig. 1. 

2 One-third subharmonic resonance and Rossby 
waves 

Figure 1 shows (as dashed line) the instrumental monthly global temperature 
data (land and ocean combined) for the period 1880-2009, as presented at 
the NOAA site http://lwf.ncdc.noaa.gov/oa/climate/ research/anomalies/ 
index.html (see also Ref. p3]). The solid curve (trend) in the figure corre- 
sponds to a wavelet (symmlet) regression of the data (cf. Refs. [T5].|16]). 
Figure 2 shows corresponding detrended fluctuations, which produce a sta- 
tistically stationary set of data. We use a nonlinear nonparametric wavelet 
regression with a relatively small number of wavelet coefficients represent- 
ing the underlying regression function. A threshold method is used to keep 
or kill the wavelet coefficients. In this case, in particular, the Universal 
(VisuShrink) thresholding rule with a soft thresholding function was used. 
Figure 5 shows a spectrum of the wavelet regression detrended data calcu- 
lated using the maximum entropy method (because it provides an optimal 
spectral resolution even for small data sets). One can see in this figure a 
small peak corresponding to a one-year period and a huge well defined peak 
corresponding to a three-years period. 

This is the one-third subharmonic resonance with the driving frequency 
oj corresponding to the annual NS-asymmetry of the solar forcing (the huge 
peak in Fig. 5 corresponds to the first term in the right-hand side of the Eq. 
(9))- 
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Figure 3: The reconstructed air temperature data (dashed line) for the pe- 
riod 0-340 kyr. The solid curve (trend) corresponds to a wavelet (symmlet) 
regression of the data. 



■ Wavelet regression detrended 




50000 100000 150000 200000 250000 300000 350000 

Years BP 



Figure 4: The wavelet regression detrended fluctuations from the data shown 
in Fig. 3. 
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Figure 5: Spectrum of the wavelet regression detrended fluctuations of the 
monthly global temperature anomaly (land and ocean combined). 
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Figure 6: The monthly global Sea Surface Temperature Anomalies (dashed 
line) for the period 1850-2008yy. The solid curve (trend) corresponds to a 
wavelet (symmlet) regression of the data. 
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Figure 7: The same as in Fig. 5 but for the wavelet regression detrended 
fluctuations of the global Sea Surface Temperature Anomalies. 
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Figure 8: Defect of autocorrelation function versus r (in In- In scales) for the 
wavelet regression detrended fluctuations of the global Sea Surface Temper- 
ature Anomalies. The straight line is drawn in order to indicate scaling Eq. 
(10). 



9 




10" 4 10" 3 10" 2 10" 1 



f [d 1 ] 

Figure 9: Spectrum of sea surface height fluctuations [18] (TOPEX/ Posei- 
don and ERS-1/2 altimeter measurements) in the logarithmic scales. The 
profound peak corresponds to the annual cycle. The dashed straight line is 
drawn in order to indicate correspondence to the scaling Eq. (10). 
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Figure 10: The monthly Global-SST ENSO index (dashed line) for the period 
1850-2008yy (the index is in hundredths of a degree Celsius). The solid curve 
(trend) corresponds to a wavelet (symmlet) regression of the data. 
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Figure 11: The same as in Fig. 5 but for the wavelet regression detrended 
fluctuations of the Global-SST ENSO index. 

The fluctuations of oceanic temperature cause certain variations of the sea 
surface height. These variations are intermixed with the sea surface height 
variations caused by the oceanic planetary Rossby waves. The oceanic plane- 
tary Rossby waves play an important role in the response of the global ocean 
to the forcing (see, for instance, Refs. [H],[l8]) and they are of fundamen- 
tal importance to ocean circulation on a wide range of time scales (it was 
also suggested that the Rossby waves play a crucial role in the initiation 
and termination of the El Nino phenomenon, see also below). Therefore, 
they present a favorable physical background for the global subharmonic res- 
onance. It should be noted that the variable h in equation (7) represents 
western boundary Rossby wave reflection of its counterpart in the equation 
(6). That provides a negative feedback to the T tendency in the eastern 
equatorial ocean via a non-linear thermocline displacement: h + bh 3 (its ver- 
tical gradient separating near-surface and deep-ocean mixed layers [3], [I]). 
The solar periodic forcing plays analogous role in this equation. 

For that reason it is interesting to look also separately at global sea surface 
temperature anomalies. These data for time range 1850-2008yy are shown in 
Fig. 6 (the monthly data are available at |http: / /jisao.washington.edu/data/global| 
sstanomts/, see also Ref. [IH]). The solid curve (trend) in the figure corre- 
sponds to the wavelet (symmlet) regression of the data. Figure 7 shows 
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Figure 12: The same as in Fig. 8 but for the wavelet regression detrended 
fluctuations of the Global-SST ENSO index. 

a spectrum of the wavelet regression detrended data calculated using the 
maximum entropy method. The spectrum seems to be very similar to the 
spectrum presented in Fig. 5. Since the high frequency part of the spec- 
trum is corrupted by strong fluctuations (the Nyquist frequency equals 0.5 
[mon -1 ]), it is interesting to look at corresponding autocorrelation function 
C(t) in order to understand what happens on the monthly scales. It should 
be noted that scaling of defect of the autocorrelation function can be related 
to scaling of corresponding spectrum: 

1 - C(r) ~ r a E(f) ~ r (1+a) (10) 

Figure 8 shows the defect of the autocorrelation function in ln-ln scales in 
order to estimate the scaling exponent a ~ 0.6±0.04 (the straight line in this 
figure indicates the scaling Eq. (10)). The existence of oceanic Rossby waves 
was confirmed rather recently by NASA/CNES TOPEX/Poseidon satellite 
altimetry measurements. Corresponding to these measurements spectrum 
of the sea surface height fluctuations, calculated in Ref. [18J with daily 
resolution (see also [20]), is shown in Figure 9. The dashed straight line in 
this figure is drawn in order to indicate correspondence to the scaling Eq. 
(10): 1 + a ~ 1.6. 
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Figure 13: Spectrum of the wavelet regression detrended fluctuations shown 
in Fig. 4. 

The Rossby waves (together with Kelvin waves) and a strong atmosphere- 
ocean feedback provide physical background for the El Nino phenomenon 
(see, for instance, Refs. [3], [I], [21] and references therein). Figure 11 shows 
spectrum for the wavelet detrended fluctuations of the so-called Global-SST 
ENSO index (Fig. 10), which captures the low- frequency part of the El Nino 
phenomenon (the monthly data are available at |http:/ /jisao. Washington .edu/| 
data/globalsstenso/). The annual forcing can come from the oceanic Rossby 
waves (cf Fig. 9). To support this relationship we show in figure 12 defect of 
autocorrelation function calculated using the wavelet detrended fluctuations 
from Fig. 10. The In- In scales have been used in Fig. 12 in order to estimate 
the scaling exponent a ~ 0.6 ± 0.03 (the straight line in this figure indicates 
the scaling Eq. (10), cf. Figs. 8 and 9). Using these observations one 
can suggest that the El Nino phenomenon has the one-third subharmonic 
resonance as a background. 

3 Nonlinear Paleoclimate 

Recent paleoclimate reconstructions provide indications of nonlinear prop- 
erties of Earth climate at the late Pleistocene [22], [23] (the period from 0.8 
Myr to present). Long term decrease in atmospheric CO2, which could re- 
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suit in a change in the internal response of the global carbon cycle to the 
obliquity forcing, has been mentioned as one of the principal reasons for this 
phenomenon (see, for instance, [21]- [26]). At present time one can recognize 
at least two problems of the nonlinear paleoclimate, which we will address in 
present paper using recent data and speculations [27]. 

A. Reconstructed air temperature on millennial time scales are known to 
be strongly fluctuating. See, for instance figures 3 and 4. While the nature of 
the trend is widely discussed (in relation to the glaciation cycles) the nature 
of these strong fluctuations is still quite obscure. The spectral analysis of 
the wavelet regression detrended data reveals rather surprising nature of 
the strong temperature fluctuations. Namely, the detrended fluctuations of 
the reconstructed temperature are completely dominated by the one-third 
subharmonic resonance, presumably related to Earth precession effect on 
equatorial insolation. 

B. Influence of Galactic turbulent processes on the Earth climate can be 
very significant for time-scales less than 2.5 kyr. 

Figure 3 shows reconstructed air temperature data (dashed line) for the 
period 0-340 kyr as presented at the site |http://www.ncdc.noaa.gov/paleo/| 
metadata/noaaicecore-6076.html (Antarctic ice cores data, see also Ref. [H]). 
The solid curve (trend) in the figure corresponds to a wavelet (symmlet) re- 
gression of the data (cf. Ref. [IS]). Figure 4 shows corresponding detrended 
fluctuations, which produce a statistically stationary set of data. Figure 13 
shows a spectrum of the wavelet regression detrended data calculated us- 
ing the maximum entropy method. One can see in this figure a small peak 
corresponding to period ~ 5kyr and a huge well defined peak corresponding 
to period ~15kyr. We also obtained analogous results (approximately 10% 
larger) from the "Vostok" ice core data for period 0-420kyr (for the data 
description see Refs. [28]. [29]). 

Origin of the periodic energy input with the period ~ 5kyr can be re- 
lated to dynamics of the energy that the intertropical regions receive from 
the Sun (equatorial insolation). Indeed, it is found in Ref. [30] that a clear 
and significant 5kyr period is present in this dynamic over last 1 Myr. The 
amplitude of the 5kyr cycle in the insolation decreases rapidly when getting 
away from the equator. Using the fact that double insolation maximum and 
minimum arise in the tropical regions in the course of one year, the authors 
of the Ref. [30] speculated that this period in seasonal amplitude of equa- 
torial insolation is determined by fourth harmonic of the Earth precession 
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Figure 14: A small-time-scales part of the autocorrelation function defect of 
the wavelet regression detrended fluctuations from the data shown in Fig. 4. 
The straight line indicates the Kolmogorov's '2/3' power law for the structure 
function. 

cycle. It should be noted, that the idea of a significant role of tropics in gen- 
erating long-term climatic variations is rather a new one (see Ref. [30] for 
relevant references). In Ref. [31], for instance, the authors speculated that 
the high frequency climate variability (in the millennial time scales) could be 
related to high sensitivity of the tropics to summer time insolation. Then, the 
oceanic advective transport could transmit an amplified response of tropical 
precipitation and temperature to high latitudes. Physical mechanism of this 
amplification is still not clear and the above discussed one-third subharmonic 
resonance can be a plausible possibility (in this respect it is significant, that 
we used the Antarctic data). 

4 Galactic turbulence and the temperature 
fluctuations 

Since the high frequency part of the spectrum is corrupted by strong fluctu- 
ations (the Nyquist frequency equals 0.5 [(500?/) _1 ]), it is interesting to look 
at corresponding autocorrelation function C{r) and at structure functions 
S p {t) (of different orders p) in order to understand what happens on the mil- 
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lennial time scales. The correlation function defect 1 — C(r) is proportional 
to the second order structure function S 2 (t). Therefore, we can compare 
results obtained by these different tools. First let us look at autocorrelation 
function C(r). Figure 14 shows a relatively small-times part of the correla- 
tion function defect. The ln-ln scales have been used in this figure in order to 
show a power law (the straight line) for the second order structure function: 
S 2 (r) = (\x(t + r) - x(t)| 2 ) : 

l-C(r) ocS 2 (t) ocr 2/3 (11) 

This power law: '2/3', for structure function (by virtue of the Taylor hypoth- 
esis transforming the time scaling into the space one [32], [33]) is known for 
fully developed turbulence as Kolmogorov's power law. 

Although, the scaling interval is short, the value of the exponent is rather 
intriguing. This exponent is well known in the theory of fluid (plasma) tur- 
bulence and corresponds to so-called Kolmogorov's cascade process. This 
process is very universal for turbulent fluids and plasmas [34j,[35j. In spite of 
the fact that magnetic field is presumably important for interstellar turbu- 
lence the Kolmogorov description can be still theoretically acceptable even 
in this area [3S]-[SH]- Moreover, the Kolmogorov-type spectra were observed 
on the scales up to kpc. In order to support the Kolmogorov turbulence as a 
background of the wavelet regression detrended temperature modulation we 
calculated also structure functions S p (r) = (\x(t + r) — x(t)\ p ) with different 
orders p. In the classic Kolmogorov turbulence (at very large values of the 
Reynolds number [3"2"].[39j) 

S p oc r<* C P ^ | (12) 

(at least for p < 3). Figure 15 shows a small-time-scales part of the struc- 
ture functions S p with p = 0.2, 0.5, 0.7, 1, 2, 3 for the wavelet regression 
detrended fluctuations from the data shown in Fig. 4. The straight lines 
are drawn in order to indicate scaling in the ln-ln scales. Figure 16 shows 
as circles the scaling exponent ( p against p for the scaling shown in Fig. 15 
(the exponents were calculated using slopes of the straight lines - best fit, in 
Fig. 15). The bars show the statistical errors. The straight line in Fig. 16 
corresponds to the strictly Kolmogorov turbulence with ( p = p/3 Eq. (12). 
One can see good agreement with the Kolmogorov turbulence modulation. 

For turbulent processes on Earth and in Heliosphere the Kolmogorov's 
scaling with such large time-scales certainly cannot exist. Therefore, one 
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Figure 15: A small-time-scales part of the structure functions S p (p=2,3) for 
the wavelet regression detrended fluctuations from the data shown in Fig. 
4. The insert shows the structure functions for p = 0.2, 0.5, 0.7, 1. The 
straight line indicates scaling in the ln-ln scales. 

should think about a Galactic origin of Kolmogorov turbulence (or turbulence- 
like processes jlO]) with such large time-scales. Let us recall that diameter 
of the Galaxy is approximately 100,000 light years. This is not surprising if 
we recall possible role of the galactic cosmic rays for Earth climate. Galactic 
cosmic ray intensity at the Earth's orbit is modulated by galactic turbu- 
lence [33]. On the other hand, the galactic cosmic rays can determine the 
amount of cloud cover (a very significant climate factor) on global scales 
through the massive aerosols formation (see, for instance, [H|-[50j). Thus, 
the galactic turbulence can modulate the global temperature fluctuations by 
the Kolmogorov scaling properties on the millennial time scales. If one knows 
the characteristic velocity scale v for the Taylor hypothesis one can estimate 
outer space-scale of corresponding galactic turbulence as L > 2500?/ x v. 
However, it is not clear what estimate we should take for the v. For instance, 
one could try velocity of the solar system relative to the cosmic microwave 
background (CMB) rest frame: v ~ 370 km/sec. In this case one obtains 
L ~ lpc. It should be noted that in recent paper it is suggested that 
the typical outer scale for spiral arms can be as small as lpc and in interarm 
regions the outer scale can be larger than lOOpc. Since the solar system and 
Earth are at present time within the Orion Arm this suggestion is in agree- 
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Figure 16: The scaling exponent ( p against p for the scaling shown in Fig. 
15 (circles). The straight line is drawn in order to indicate the Kolmogorov's 
turbulence with ( p = p/3 Eq. (12). 

ment with the above estimate. Although the suggestion of the Ref. jl9] is 
still under active discussion the paleoclimate consequences of this suggestion 
can be very interesting and we will discuss one of them here. Namely, when 
orbiting the Galactic center the solar system and Earth are in the interarm 
regions the reverse of the Taylor hypothesis provides us with the outer time- 
scale ~ 2500 x 100 years. This time scale is larger then any known glaciation 
period (which are determined by the periods related to orbiting Earth around 
Sun, see for instance [7]). Strong turbulent fluctuations of the cosmic rays 
flux on such large time-scales should prevent to the glaciation cycles to occur 
when the solar system is in the interarm regions. The Earth deglacitaion 
related to the interarm regions was suggested in Refs. [I5],[l6] and explained 
by difference in intensity of the cosmic ray flux in the spiral arms and in the 
interarm regions. It is difficult to estimate at present time which of the two 
mechanisms is more efficient. In any way they are working to the same end 
and both are open to discussion. 
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5 Subharmonic resonance in solar activity 



The solar activity is chaotic but has a well-defined mean period of about 11 
years. The 11-year cycle is well known for more than a century and a half. 
Despite this, nature of the 11-year cycle is still a subject of vigorous inves- 
tigations. The most popular point of view is a 'dynamo-wave' mechanism. 
It is assumed that a magnetic dynamo, generated by the solar differential 
rotation and the helicity of turbulent convective flows, produces this propa- 
gating wave. Helicity, through the a-effect, plays crucial role in this mecha- 
nism. Recently, observational data of physical quantities associated with the 
a-effect became available and a considerable progress in this direction was 
achieved (see, for instance, Refs. [SI] , [52] , [53] ) • The a-effect has two contrib- 
utors: one related to helicity of convective vortices and another related to 
the helicity of magnetic field. Intrinsically non-linear character of the prob- 
lem makes it especially difficult for theoretical investigation. The nonlinear 
solar dynamo has to be saturated in order to get a quasi-stationary wave. 
The magnetic part of the a-effect can play a crucial role in such saturation, 
while certain modification of the turbulent diffusivity and other transport 
coefficients is unavoidable at this process (see Refs. [51] , [52] , [53] and sec- 
tion 'Chaotic dynamo'). Results of the above mentioned simulations show 
that the dynamo model leads to a steadily oscillating magnetic configura- 
tion. The cyclic behavior is typical for moderate dynamo action whereas 
for the stronger dynamo action a chaotic behavior is usually observed in the 
dynamo simulations. It is believed that the dynamo mechanism is mainly 
operational deep inside the convective zone (or even in the overshoot layer). 
However, in the upper layers of the convection zone ('at surface') a strong 
hydromagnetic activity has been also observed. The presence of large-scale 
meandering flow fields (like jet streams), banded zonal flows and evolving 
meridional circulations produces a very complex picture [51] . In this situ- 
ation one can expect that more than one non-linear mechanisms can be at 
background of the observed 11-year solar cycle. Namely, it could be inferred 
from the data analysis presented below that a dynamo mechanism, directly 
affected by the solar rotation, may generate the basic chaotic oscillations 
(with a fundamental period different from the 11-year cycle), while another 
non- linear mechanism amplifies these oscillations to the observed 11-year 
chaotic oscillations. In this two-stage picture the mirror asymmetry of the 
solar magnetic field can be still a key driver of the 11-year activity cycle 
through an additional non- linear amplifying mechanism. 
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Figure 17: The monthly sunspot number (dashed line) for the period 1749- 
2009 years. The solid curve (trend) corresponds to a wavelet (symmlet) 
regression of the data. 
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Figure 18: The wavelet regression detrended fluctuations from the data 
shown in Fig. 17. 
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Most of the regression methods are linear in responses and statistical 
analyses of the experimental sunspot data was dominated by linear stochastic 
methods, while it was recently rigorously shown in Ref. (55] that a nonlinear 
dynamical mechanism (presumably a driven nonlinear oscillator, see also Ref. 
|56j) determines the sunspot cycle. Figure 17 shows the monthly sunspot 
number (dashed line) for the period 1749-2009 years (the data are available 
at |http: / / sidc.oma.be/ sunspot-data/ ) . The solid curve (trend) corresponds 
to a wavelet (symmlet) regression of the data (cf. Refs. [IS]). Figure 18 
shows corresponding detrended fluctuations, which produce a statistically 
stationary set of data. The wavelet detrended data set is not statistically 
stationary, however, for the extended 1611-2009 period, presumably due to 
the Maunder minimum in the 17th century (see Refs. [57], [58]). Figure 
19 shows a spectrum of the wavelet regression detrended data calculated 
using the maximum entropy method (because it provides an optimal spectral 
resolution even for small data sets). In Fig. 19 one can see a well defined 
peak corresponding to period ~ 3.7 years. 

The wavelet regression method detrends the data from the approximately 
11-years period (cf. Fig. 17). Therefore, it is plausible that the one-third 
subharmonic resonance [8| can be considered as a background for the 11- 
years solar cycle: 11/3.7 ~ 3. Indeed, it is known [59J that interaction of the 
Alfven waves (generated in a highly magnetized plasma by a cavity's moving 
boundaries) with slow magnetosonic waves can be described using Duffing 
oscillators (see also Refs. [55], [60] and below). 

6 Chaotic dynamo 

In order to understand appearance of the 3. 7- years period let us represent the 
spectrum shown in Fig. 19 in semi-logarithmical scales: figure 20. In these 
scales an exponential behavior corresponds to a straight line. It is known, 
that both stochastic and deterministic processes can result in the broad-band 
part of the spectrum, but the decay in the spectral power is different for the 
two cases. The exponential decay indicates that the broad-band spectrum for 
these data arises from a deterministic rather than a stochastic process. For a 
wide class of deterministic systems a broad-band spectrum with exponential 
decay is a generic feature of their chaotic solutions Refs. [61]- [M]- A wavy 
exponential decay (see Fig. 20) is a characteristic of a chaotic behavior 
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Figure 19: Spectrum of the wavelet regression detrended fluctuations shown 
in Fig. 18. 
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Figure 20: The same as in Fig. 19 but in semi-logarithmical scales. The 
dashed straight line indicates an exponential decay. 
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generated by time-delay differential equations [22] • A classic example of 
time-delay differential equation with chaotic solutions is the Mackey-Glass 
equation: 

du{t) 0.2 -u(t-r) 

— ; — — tz 7 rr^r — 0.1 ■ u[t) 13 

dt (l + w(t-r) 10 ) w y J 

Figure 21 shows spectrum of a solution of this equation for the time-delay 
r = 30. The dashed straight line indicates an exponential decay (cf. Fig. 
20). 

In the Parker dynamo, which was generalized in the Refs. [51] , [52] , [53] , a 
time-delay in the back influence of the magnetic field on the a-effect [65j,[66j 



(14) 

(where B is the azimuthal component of the magnetic field) can significantly 
change the evolution of the magnetic field even for a small time delay. In 
particularly, this non-linear delay can result in appearance of processes with 
periods much longer than the fundamental period through a parametric res- 
onance [66J. 

In the dynamo models that have physically distinct source layers the 
finite time is required in order to transport magnetic flux from one layer 
to another (a time-delay involved in the a-quenching mechanism due to the 
Lorentz feedback). Even the Duffing equations for B can be obtained for the 
dynamo models using the delay idea in this case (cf. Eq. (14) in the Ref. 

my- 

B + tu 2 B + 1 B + /3B(t-r)f(B(t-r)) = (15) 

where f(B(t — r)) is a quenching factor, which can be approximated by a 
nonlinear function. In particular, in the Ref. [67] (see also Ref. [68]) this 
function has been approximated as: 

f(B) = 1(1 + erf (B 2 - B 2 mm )(l - erf (B 2 - B 2 mm ) (16) 

and for small (B 2 — B 2 min ) we obtain a nonlinear delayed Duffing equation 

B + Vl 2 B + -fB + (3'B 3 (t-r) = (17) 

where f2g,7 and 0' are certain constants. The subharmonic resonances and 
chaotic regimes are also known for the delay Duffing equations with a periodic 
forcing. 
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Figure 21: Spectrum of the Mackey-Glass chaotic time-series. The dashed 
straight line indicates an exponential decay. The insert shows a sketch of 
corresponding complex-time plane. 

It is also significant for those dynamo models that have spatially segre- 
gated source regions for the poloidal and toroidal magnetic field components 
(such as, for instance, the Babcock- Leighton dynamo mechanism (68]). In 
the global dynamo models that include meridional circulation the time delay 
related to the circulation should be comparable to global rotation period (see 
below). 

Nature of the exponential decay of the power spectra of the chaotic sys- 
tems (Figs. 20 and 21) is still an unsolved mathematical problem. A progress 
in solution of this problem has been achieved by the use of the analytical con- 
tinuation of the equations in the complex domain (see, for instance, [64J ) . In 
this approach the exponential decay of chaotic spectrum is related to a singu- 
larity in the plane of complex time, which lies nearest to the real axis (see the 
insert in Fig. 21). Distance between this singularity and the real axis deter- 
mines the rate of the exponential decay. For many interesting cases chaotic 
solutions are analytic in a finite strip around the real time axis. This takes 
place, for instance for attractors bounded in the real domain (the Lorentz 
attractor, for instance). In this case the radius of convergence of the Taylor 
series is also bounded (uniformly) at any real time. 

Let us consider, for simplicity, solution u(t) with simple poles only, and 
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to define the Fourier transform as follows 



■T e /2 



dt e- i2nft u(t) 



T e /2 



(18) 



Then using the theorem of residues 



u(f) = *(2tt) 



' 2 '£R j exp(i2Trfx j -\2Trfy j \) 



(19) 



3 



where Rj are the poles residue and Xj + iyj are their location in the relevant 
half plane, one obtains asymptotic behavior of the spectrum E(f) = \u(f)\ 2 
at large / 



where y m in is the imaginary part of the location of the pole which lies nearest 
to the real axis. In the case of symmetric analytic strip with a width A = 



(cf. the insert in Fig. 21). 

The chaotic spectrum provides two different characteristic time-scales for 
the chaotic system: a period corresponding to fundamental frequency of the 
system, ff un d, and a period corresponding to the exponential decay rate, 2ttA 
(cf. Eq. (21)). The fundamental period can be estimated using position of 
the low-frequency peak (cf. Figs. 20 and 21), while the exponential decay rate 
period 2nA can be estimated using the slope of the straight line of the broad- 
band part of the spectrum in the semi-logarithmic representation. In the case 
of the global solar dynamo the width of the analytic strip A can be theoret- 
ically estimated using the Carrington solar rotation period: A ~ T c ~ 25.38 
days. This period roughly corresponds to the solar rotation at a latitude of 
26 deg, which is consistent with the typical latitude of sunspots (cf. Fig. 20). 

Additionally to the exponential spectrum (Fig. 20), let us check the 
chaotic character of the wavelet regression detrended fluctuations calculating 
the largest Lyapunov exponent: \ ma x- A strong indicator for the presence 
of chaos in the examined time series is condition X ma x > 0. If this is the 
case, then we have so-called exponential instability. Namely, two arbitrary 
close trajectories of the system will diverge apart exponentially, that is the 
hallmark of chaos. To calculate X max we used a direct algorithm developed by 
Wolf et al. [69]. Figure 22 shows the pertaining average maximal Lyapunov 




(20) 



E(f) ~ exp(-27rA /) 



(21) 
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Figure 22: The pertaining average maximal Lyapunov exponent at the per- 
taining time, calculated for the same data as those used for calculation of the 
spectrum (Figs. 19 and 20). The dashed straight line indicates convergence 
to a positive value. 

exponent at the pertaining time, calculated for the data set shown in Fig. 
18. The largest Lyapunov exponent converges very well to a positive value 
\ mnr ^ 0.286 morT 1 > 0. 



7 Flip-flop phenomenon 

It should be noted that the same period ~ 3.7 years was recently found for 
the so-called flip-flop phenomenon of the active longitudes in solar activ- 
ity [TO], [71]. Sunspots are tend to pop up preferably in certain latitudinal 
domains and move toward the equator due to the 11-year cycle. Recently, 
strong indications of non-uniform longitudinal distribution of sunspots (ac- 
tive longitudes) was reported and analyzed in a dynamic frame related to the 
mean latitude of sunspot formation, in which the active longitudes persist 
for the last eleven solar 11-years cycles (see Refs. [70], [71] and references 
therein). At any given time, one of the two active longitudes (approximately 
180° apart) exhibits a stronger activity - dominance. Observed alternation of 
the active longitudes dominance in 3.7 years on average was called as flip-flop 
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phenomenon [7U] . It seems rather plausible that the observed flip-flop period 
and the fundamental period of the wavelet regression detrended fluctuations 
of solar activity (Fig. 20) have the same origin. In this vein, the observation 
[70], [72] that the period of the flip-flop phenomenon follows to variations of 
the real length of the sunspot cycle (which has the 11-years period on av- 
erage only) supports the idea of the one-third subharmonic resonance as a 
background of the 11-years cycle of solar activity. 

Another relevant example of the 3.7 years period appearance is the inter- 
planetary magnetic field polarity variations [73] • Close value (~ 3.5 years) 
of a period of the geomagnetic aa index was reported in Ref. [71] . 



8 Multidecadal coherence 

The question: Whether there is a coherence between solar activity and global 
temperature dynamics, was widely discussed in relation to the global warm- 
ing. Naturally, in relation to the global warming just trends (solid curves in 
the Figs. 1 and 17) of the corresponding dynamical processes were studied. 
Now we can compare the wavelet regression detrended components of these 
processes (Figs. 2 and 18). Since the dominating periods of these detrended 
dynamical processes are different (3 and 3.7 years) we will look at a low- 
frequency domain (with frequency / < 0.3?/ -1 ). The cross spectrum Ei^(f) 
of two processes x±(t) and X2(t) is defined by the Fourier transformation of 
the cross-correlation function normalized by the product of square root of 
the univariate power spectra E\(f) and -E^Cf): 



the bracket (...) denotes the expectation value. The cross spectrum can be 
decomposed into the phase spectrum (fixpif) an d the coherency Ci )2 (/): 



Because of the normalization of the cross spectrum the coherency is ranging 
from Ci )2 (/) = 0, i.e. no linear relationship between Xi(t) and x 2 (t) at /, to 
Ci^if) = 1) i-e. perfect linear relationship. 




J2 T (xi{t)x 2 {t - t)) exp(-i27r/r) 



(22) 




EiM) = Ci, 2 (/)e-^' 2(/) 



(23) 
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Figure 23: Coherency between the wavelet regression detrended solar activity 
and global temperature in a low-frequency domain. 

Figure 23 shows the coherency between the wavelet detrended solar activ- 
ity and global temperature (for the period 1880-2009yy) in the low-frequency 
domain, computed using the fast Fourier transformation. One can see a very 
high coherency for the time periods larger than 20 years. Let us recall that 
the 22y period corresponds to the Suns magnetic poles polarity switching 
and can be a base for a multidecadal chaotic coherence between solar activ- 
ity and global climate (cf. Refs. [7]). It is knowm that there is no clear 
correlation between the data dominated by the trends (shown in the Figs. 1 
and 17). Therefore, one can conclude that just the trends are presumably 
not correlated in this case, whereas the large-scale fluctuations (wavelet re- 
gression detrended: Figs. 2 and 18) exhibit a very high coherency for the 
time periods larger than 20 years. Thus, just the nonlinear nature of the 
Sun-climate interaction is the main cause of the well known difficulties in 
estalishing and investigation of this interaction. 

9 Conclusion 

In the considered complex natural systems the fluctuations are often of the 
same order as the trend itself. This phenomenon is related to a nonlinear 
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character of the fluctuations. Has this nonlinear phenomenon an universal 
nature? The main three points of the above consideration are: 

a) the known (crude) models can be reduced to the Duffing oscillator, 

b) the one-third subharmonic resonance is a generic property of this os- 
cillator, 

c) the data indicates crucial role of this resonance in generation of the 
large-scale fluctuation. 

The next generation of the models of global climate should take into ac- 
count the highly nonlinear nature of the Sun-climate interaction, which is not 
realized directly through the trends but through the large-scale fluctations 
(both for monthly and for multidecadal time scales). 

At this stage of the data accumulation and analysis limitations of this 
conclusions are unknown. However, we believe that the above considered 
results could attract attention of the nonlinear scientists to the natural sys- 
tems. Because, despite of their overwhelming complexity, these systems are 
governed by the simple physical laws. 
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